Flux creep in type-II superconductors: self-organized criticality approach 
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We consider the current density distribution function of a flux creep regime in type-II super- 
conductors by mapping the flux creep process to the dynamics of a model with a self-organized 
criticality. We use an extremal Robin Hood type model which evolves to Been's type critical state 
to treat magnetic flux penetration into a superconductor and derive an analog of the current-voltage 
characteristics in the flux creep region. 

PACS numbers: 74.60. Ec, 74.60. Ge 



I. INTRODUCTION 

In the presence of currents vortex structure in a type-II 
superconductor is subjected to Lorentz force. The value 
of this force per unit length of a vortex is = j x (fro, 
where j is the current density, B is the magnetic field, 
0o = <AjB/_B, and 4>q is the flux quantum. The Lorentz 
force acting on a unit volume of a vortex matter is there- 
fore given by Ft, — nfj, = j x B, where n(r) = B{r)/4>o 
is the density of the vortices^. 

Consider as an illustration a superconducting slab par- 
allel to the y, z-plane and assume that a certain mag- 
netic field B a is applied along the z-axis as shown in 
Fig. ^ I n this case we have j = j y and B = B z 
related by Maxwell's equation dB/dx = ~Hoj- The 
Lorentz force Ft, = jB is therefore proportional to the 
density of vortices n and its gradient dn/dx as follows 
from Fl — j B oc B dB/dx oc n dn/dx. 

In high-current density superconductors defects of the 
crystalline structure "pin" the vortices. This pinning 
leads to creation of various vortex configurations with 
flux "hills" and flux "valleys" between these hills. In- 
deed, single vortices or bundles of vortices redistribute 
spatially if the Lorentz force Ft, oc j cx dn/dx overcomes 
the pinning force F p i n . This means that n dn/dx oc -Fpinj 
i.e., a steady vortex structure consists of a composition of 
fragments with a certain slope dn/dx ^ (the flux hills). 
On a macroscopic level we have Ft, oc j and therefore the 
local depinning of vortices happens when the local cur- 
rent density j exceeds a certain critical value j c oc F p j n . 
It was pointed out by de Gennes that the flux hills with 
slopes statistically fixed by the critical current density j c 
are very much alike to the sand pilesi 

The concept of approaching to flux statics and dynam- 
ics in superconductors with high density of pinning cen- 
ters was first suggested by Been 3 .. The famous Bean's 
model assumes that the current density j is equal to the 
critical current density j c everywhere where there are cur- 
rents. Initially Been considered a magnetic field indepen- 
dent j c . As a result the spatial increase or decrease of the 
field inside the sample is linear. In particular, in a slab 
with thickness 2d (see Fig. ^| the assumption j c — Const 
leads to B = B a — HqJ c \x =p d|, where B a is the mag- 



FIG. 1: Series of Been's critical states in a slab parallel to 
the j/z-plane. A zero-field cooled sample was subjected to a 
monotonically increased field parallel to the z-axis. The slope 
of B(x) is proportional to the critical current density j c . 



netic field at the sample surface. In general, the field 
dependence of j c has to be taken into account in order 
to give a better description of the whole wealth of the 
experimental data available 4,5 . 

In this paper we treat the low-temperature flux creep 
in type-II superconductors in the framework of the 
self-organized criticality approach 6 : 7 : 8 : 9 : 1 ^ 11 ! 1 ^ 13 : 1 ^ 15 ! 16 . 
The distribution function of the current density j is 
considered in detail for an extremal Robin Hood type 
process 9 -. We show that this process results in a self- 
organized Been's type critical state with a complex dy- 
namics, which can be mapped to the low-temperature 
flux creep. This mapping is used to obtain an analog 
of a current-voltage characteristics of a type-II supercon- 
ductors in the flux creep regime. 

The paper is organized as follows. In Sec. II we discuss 
the dynamics of the low-temperature flux creep regime 
in type-II superconductors. The self-organized critical- 
ity of an extremal Robin Hood type process is treated 
in Sec. III. In Sec. IV we introduce the low-temperature 
flux creep model and derive the distribution function of 
the current density j. The magnetic flux penetration into 
a superconductor and an analog of current-voltage char- 
acteristics of a type-II superconductor in the flux creep 
regime are treated in Sec. V. Sec. VI summarizes the ob- 
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tained results. 



II. FLUX CREEP IN SUPERCONDUCTORS 

There are few mechanisms resulting in depinning of 
vortices for currents with a density less than the critical 
value. In particular, both thermally activated depinning 
and quantum tunnelling^ii result in vortices and bun- 
dles of vortices jumping from one group of pinning cen- 
ters to another. This type of vortex motion is called flux 
creepi&. The probability of depinning of a bundle of vor- 
tices strongly depends on the current density and tends 
to unity when the current density tends to the critical 
current density, i.e., j — > j c . The high level of correla- 
tions between vortices in a vortex matter leads to a very 
complicated collective behavior of vortices especially in 
thin superconducting films where the stray fields result 
in nonlocality of the problem^. Activation of one vor- 
tex can ignite a local avalanche-type motion of bundles 
of correlated vortices, i.e., in the flux creep region vortex 
matter is a system with avalanche driven dynamics. Such 
systems with avalanche driven dynamics are the subject 
of the modern theory of self-organized criticality which 
revealed a variety of power-law distributions for the spa- 
tial size and temporal duration of the avalanches^. 

Motion of vortices is accompanied by the field and cur- 
rent energies dissipation. The heat release during rear- 
rangements (flux avalanches) in a vortex matter can heat 
up the whole sample or a part of it to a temperature 
higher than the critical temperature T c . Under certain 
conditions an avalanche of even a small group of vortices 
can trigger a run-away magneto-thermal instability caus- 
ing the superconducting-to-normal transition—. 

The original Bean's model means that in the critical 
state the dependence of j on E is a highly nonlinear step- 
wise function 

._. /0, if £ = 0, m 
J_Jc \l, if E^O, [> 

where j c = j c e and e = E/E is a unit vector parallel to 
the electric field E induced by flux motion. 

It is well established now that in the narrow vicinity of 
the critical current (\j — j c \ <C j c ), i-e., in the flux creep 
regime the dependence of j on E is a very steep function 
given by the power law 



where n > 1 is a parameter and the field Eq defines 
the critical current density. It is common to define j c as 
the current density at Eq = 10~ 6 Vcm -1 . It is worth 
mentioning that for n>lwe can rewrite Eq. in the 
following logarithmic form 

j=j c + J ^lnf^), (3) 

n \E a J 



where the omitted terms are of order \/n? <C 1. 

The dependence of j on E given by Eq. © was first 
derived in the framework of the Anderson-Kim model? 
considering the thermally activated uncorrelated hopping 
of bundles of vortices. The vortex-glass^I and collective 
creep22i2i models result in more sophisticated dependen- 
cies of j on E. However, these dependencies coincide 
with the one given by Eq. J2J if j —j c *C j c - The recently 
developed self-organized criticality approach to the crit- 
ical stat o 10 ' 11 also leads to Eq. © if j — j c <C j c - In 
the interval j — j c <C j c the power law © is in good 
agreement with numerous experimental data™. 

A detailed study of the dynamics of vortices in the flux 
creep regime was performed by Field and coworkers™ In 
their experiments the magnetic field outside a tubular su- 
perconducting sample is ramped slowly, driving the flux 
into the tubes outer wall. After the flux front reaches 
the inner wall of a tube it spills out into the tube's in- 
terior. The entrance of flux into the tube's interior was 
detected in real time. This experiment distinguished be- 
tween flux leaving the superconductor in discrete bundles 
or avalanches. It was shown that the probability D(s) of 
an avalanche containing s vortices is a power-law depen- 
dence extending over 1.6 decades. 

Been's critical state dynamics is typical to some other 
spatially extended dynamical systems. High number of 
degrees of freedom in these systems introduces the prob- 
lem of treating the effect of coupling between the in- 
dividual degrees of freedom. In many cases even very 
complicated systems "self-organize" so that their descrip- 
tion reduces to a small number of collective degrees of 
freedom— V 

In some dynamical systems the individual degrees of 
freedom keep each other in a stable balance, which cannot 
be described as a "perturbation" of some decouple state 
as well as in terms of a small amount of collective degrees 
of freedom. This type of the self-organized systems has 
to be quite robust. In the opposite case these systems 
would not be able to evolve to a stable balanced "critical" 
state§i The sand piles and flux hills in superconductors 
demonstrate many features typical for the self-organized 
critical state. 

Several different models were suggested to illustrate 
and study dynamical systems with extended spatial de- 
grees of freedom and many metastable states. These sys- 
tems evolve into a self-organized critical state without a 
detailed specification of the initial conditions, i.e., the 
critical state is an attractor of their dynamics which is 
robust with respect to variations of parameters and the 
presence of quenched randomness. 



III. SELF-ORGANIZED CRITICALITY OF 
EXTREMAL PROCESSES 

We consider now a specific subclass of extremal pro- 
cesses demonstrating the self-organized criticality. In ex- 
tremal models at each step only the sites satisfying a 
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certain extremal criterion are involved into the system 
dynamics. In particular, the problem of low-temperature 
flux creep in superconductors can be mapped to the 
Robin Hood type extremal model 9 . 



A. Robin Hood processes. 

The following picture illustrates the main idea of the 
Robin Hood type process. The famous robber leaves ev- 
eryday his forest for a nearest market and tries to obtain 
justice. He finds the richest seller, takes from him some 
random amount of goods or money and distributes it 
among neighbors. He takes nothing to himself so that 
overall amount of goods on the market is conserved. The 
market reaches a critical state after some transition time. 

Many physical systems with a dynamical variable 
matching a certain conservation relation can be treated 
in the framework of the Robin Hood model. In partic- 
ular, the low-temperature flux creep, motion of disloca- 
tions through a random array of point obstacles, domain 
walls in ferromagnets, grain boundaries in polycrystals, 
etc, can be mapped to this type of extremal processes 9 . 
We will use the Robin Hood type model to treat the self- 
criticality of the low-temperature flux creep in supercon- 
ductors in detail. 



B. Low-temperature flux creep 

If the temperature is very low than vortices in sites 
with even slightly different current density have strongly 
different probability of depinning. Therefore, the depin- 
ning of vortices will happen first in the place with the 
highest current density. 
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FIG. 2: Illustration of the dynamics rules of the Robin Hood 
type low-temperature flux creep model. 

The model suggested by Zaitse^ simulates a part of 
a superconductor with L sites using the closed boundary 
conditions. Numbers j, (0 < i < L) describe the values of 
the current density j along the x-axis (the current flows 



along the y-axis). At each simulation step a site number 
m with the maximal current density j m is found (see 
Fig. This current density is reduced by Aj chosen 
from a uniform distribution U(j) 

jm —> jm - Aj, j m ±i -> j m ±i + Aj/2. (4) 

It is worth mentioning that Zaitsev's model conserves the 
total current 

L 

I = J2jn=(j)L, (5) 

71=1 

which means that the average value of the critical current 
density (j) is also conserved. The system approaches the 
critical state starting from any initial state with a given 
value of (j). In this critical state almost all sites have cur- 
rent densities less than a certain value j c . Each site of 
the system maps an area in a superconductor containing 
a large amount of vortices but depinning is referred to a 
single vortex. Motion of this vortex leads to a rearrange- 
ment of a cluster of correlated vortices and changes the 
local current density by some value which is determined 
by the uniform distribution U(j). 

The distribution of the values of Aj which will be re- 
ferred below as A., -distribution, defines the number of 
vortices participating in an avalanche starting from a de- 
pinning of a single vortex, or, to be more precise, it is 
related to the number of vortices leaving a given site as 
a result of an avalanche. At sufficiently low temperature 
this vortex rearrangement is a very fast process compar- 
ing to the thermally activated depinning. Therefore, all 
the details concerning these avalanches, including micro- 
scopic vortex dynamics, are "hidden" in A 3 -distribution. 

C. Dynamics of extremal processes 

The fact that the flux creep dynamics conserves the 
total current / leads to a high level of correlation be- 
tween the sites. In the extremal models the extremal 
sites are located according to some problem-specific cri- 
terion and they provoke changes of the neighbor sites, 
again, according to rules specific for a given model. We 
will call these sites ignition sites. The sites drawn in the 
activity by the ignition sites will be called involved sites. 
An involved site in its order can became a ignition sites 
if it matches a certain extremal criterion. A system in 
a critical state is characterized by a critical value of the 
dynamical parameter. In particular, for the flux creep 
model this parameter is j c which is the least value of j 
in the ignition sites. That means that for the flux creep 
model all ignition sites have values greater than j c . 

The rules of dynamics given by Eq. @ are illustrated 
schematically in Fig. |3J The interval of values of current 
densities in the ignition sites we call active zone. Each 
site in the active zone, active site, at a certain moment 
becomes an ignition site. At any step only a small part 
of all sites is active. As shown in Fig. [31 the majority of 
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the sites belongs to the calm zone. We will consider the g <j> 
dynamical properties of the low-temperature creep model 03 
(@J using this terminology. 
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FIG. 3: The calm and the active zones for a one-dimensional 
flux creep model. The rules of dynamics given by Eq. Q arc 
illustrated by the dashed lines. 
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FIG. 4: The dependence of the distribution function G(j) on 
the current density j. 



IV. LOW TEMPERATURE CREEP MODEL 

The results of our numerical simulations of a system 
obeying the dynamics rules given by Eq. |0J demonstrate 
that G(j) is described by an exponent up to a certain 
critical value j c , i.e., 

G(j) = A exp(j/j e ), for j c - 1 < j < j c , (6) 

where A and j e are parameters of the distribution (see 
Fig. 0J. The distribution function G(j) has a sharp cut 
off at j = j c . In the active zone (j > j c ) the "tail" of 
G(j) decreases as 1/L, when L — > 00 (see Fig. |5j|. 

The distribution function G(j) has another cut off at 
j c — 1 if the Aj -distribution in the interval (0, 1) is cho- 
sen to be uniform. Indeed, there are only two options to 
create a site with a current density j: (a) to decrease the 
current density in the ignition site by Aj and (b) to in- 
crease the current density in one of the involved sites by 
Aj/2. The lowest value of j is obtained by subtracting 
from the minimum value of the current density in an igni- 
tion site, which is j c , the maximum value from U, which 
is 1. As a result the distribution function G(j) has a cut 
off from the left. This cut off is not universal and de- 
pends on the Aj -distribution. A continuously decreasing 
Aj-distribution eliminates both the left cut off and the 
deviation of G(j) from the exponential function in the 
vicinity of j c . The details of a specific Aj-distribution 
affect only the tail of the function G(j), i.e., its behavior 
in the interval j > j c . Therefore, we conclude that the 
values of the current density j are distributed according 
to the exponential law given by Eq. © in almost all sites 
of the system. 



A. Relation between the parameter j e and j c 

The relation between the parameter j e and the critical 
current density j c can be calculated analytically using the 



normalization and average current density conservation 
rules. Indeed, normalization of the distribution function 
G(j) given by Eq. © 

f C G(j) dj = A r exp(i/i e ) dj = 1 (7) 
J — 00 J — 00 

relates A, j c , and j e 

A=- exp(-j c /j e ). (8) 

Je 

Next, we calculate the average value of the current den- 
sity 

(i>= f C jG(j)dj (9) 

J —OO 

and obtain the dependence of A, j c , and j e on (j) 

(j) = Aj 2 e (jc/je ~ 1) exp(i c /j e ). (10) 

Combining Eqs. JSJ and i|l(J|) we find the relation be- 
tween the critical current density j c , parameter of the 
exponent j e , and average current density (j) 

3c-3e=(j). (11) 

If Aj-distribution is uniform, then we have (j) = 1/2 
and we find 

j c - Je = 1/2. (12) 

Our numerical simulations show that the exponential 
form 10 of the distribution function G(j) is universal, 
meaning that it does not depend on the spatial dimen- 
sion of the model, on the form of the A^-distribution, on 
the dynamic rules of redistribution of Aj between neigh- 
bors, and on the decision on how many sites from the 
vicinity of an ignition site are involved in the redistribu- 
tion process (for example, next nearest neighbors can be 
included into the dynamics). Therefore, we conclude that 
the exponential behavior of G(j) is typical for extremal 
models with short range interactions and a conservation 
relation for the dynamical variable. 
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Using Stirling's formula we have for a 
a = L In L — rt,; In r. 



lnT 



(15) 



The maximum of a can be found using Lagrange method 
for conditional extremum with two conditions Ijl^jl 



F = a + (a + 1)L - 0{j) L 
dF 

- — = - In m + a + /3ji = 



and we get finally 



with a and (3 defined by 



(j)L. 



(16) 



(17) 



(18) 



FIG. 5: Distribution functions G(j) for different values of the 
length of the system L in the region j > j c ["tails" of G(j)]. 



B. Origin of the exponent 

We performed numerical simulations for different ex- 
tremal models and obtained the G(j) distribution for all 
of them. In addition to the uniform distribution for Aj 
we tested also the exponential, gaussian, and power law 
distributions. These A 3 ; -distributions are more "natural" 
than U(j), since they decay gradually and have no cut- 
offs, unlike the uniform distribution. All of them lead to 
the exponential form of G(j) if two conditions are sat- 
isfied, (a) Interactions in the model are local, meaning 
that an activated site affects only sites in its surrounding 
of a finite size; (b) At each step the sum of the dynamical 
variables (ji) stays constant. 

The exponential behavior of G(j) can be obtained ana- 
lytically as the most probable distribution of the dynam- 
ical variable ji. To simplify the derivation we use discrete 
values of j dividing the whole interval of values of j into 
small intervals. In this case each site is characterized by a 
certain value ji. We denote the number of sites with the 
same value ji as rii and write the conservation relations 
in the form 



(13) 



Distribution of L sites between the intervals with j = ji 
is described by a set {rii} of numbers of sites in each 
interval. The number of states corresponding to the same 
set {rii} is given by 



r = 



L! 



n, 



(14) 



Using the continuous form of Eq. i|18|) 
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je a+ ^dj = {j)L (19) 



we find 



/?=-, e a =A=- exp(-j c /j e ), (20) 

Je Je 

i.e., the values j c , je, and (j) are related by Eq. Ijlljl and 



G(j) 



Je 



exp 



3 — Jc 

je 



1 f j (j) 

— exp ; 

eje V Je 



(21) 



The distribution function given by Eq. I|21|) contains 
only one independent parameter, j e which is character- 
izing the width of the distribution function G{j). The 
value of je is proportional to the width SAj of the Aj- 
distribution (j e oc SAj) and depends, in particular, on 
the spatial dimension of the model, on the number of the 
neighbors of the ignition site, etc. 

It is worth mentioning that the main features of the 
distribution function G(j) can be formulated in terms 
of a certain effective "temperature" associated with the 
system. Indeed, the above derivation of the distribution 
function G{j) is based on the same arguments which are 
used to derive the Gibbs distribution for a system of non- 
identical particles. We can write Eq. I|21|l in the form 
G(j) — exp(—€j/T)/j e , where the effective "energy" ej 
and "temperature" r are defined as = a (j c — j) and 
t = aj e and a is an arbitrary coefficient. Using Eq. i|ll|) 
we find for the average energy 



(ej) =a(j c - (j)) = aj e 



(22) 



The value of the parameter j e is proportional to the width 
SAj of the Aj -distribution, and therefore the effective 
temperature r tx SAj . 
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V. FLUX PENETRATION MODEL 

Extremal models can be useful for studying macro- 
scopic processes in superconductors and, in particular, 
the process of penetration of magnetic flux. In this sec- 
tion we modify our model to allow for considering this 
process. The periodic boundary conditions which we 
used above simulated an inner part of the system located 
far from the sample edges. 

Consider now a slab subjected to a magnetic field par- 
allel to its surface as shown in Fig. ^ The symmetry of 
the problem allows to perform the numerical simulations 
in one half of the sample. Assume that the applied field 
rises up with a certain rate h and the flux penetrates in- 
side the sample from its left (x = —d) edge. In this case 
the boundary condition at x = — d takes the form 

j {t + St)=j {t) + h5t, (23) 

where jo(t) is the current density at x = —d. 

The dynamics rule at the middle of the slab (x = 0) 
has the form 

3L(t + St) = j L {t)-A 3 , 
j L -i(t + St) = j L ^(t) +Aj/2, (24) 

where jh{t) is the current density at x — 0. It is worth 
mentioning that there is no current conservation at the 
sample boundaries at each simulation step. However, in a 
stationary state the current conservation holds in average 
for large time intervals. 
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FIG. 6: Starting from the same initial distribution, the cur- 
rent density approaches different asymptotic values depending 
on the ramp rate h = 4xl0" 5 , 3xl0" 5 , 2xl0~ 5 , 5xl0" 6 , 3x 
10~ 6 (from top to bottom). 



The mean time between two depinning events (5t) is in- 
versely proportional to Pd and therefore we write 



(St)(t) = Stack exp 



jm(t) 

h 



(26) 



where (fatick is the time interval corresponding to one tick 
of the "real" time clock. The last equation provides a 
synchronization mechanism between the steps in the nu- 
meric simulations and the "real" time intervals. 

Using Eqs. (|23|l and l|26|l we arrive to the synchronized 
boundary condition at x = —d in the form 



A. Introduction of the "real" time scale 

The above self-organized criticality approach to the 
process of the low-temperature flux creep formulates the 
dynamical rules of the current density dynamics and op- 
erates with the simulation steps only. These steps corre- 
spond to the sequential depinning events in the process 
of numerical simulation of the model. The "real" time 
elapsed between two successive steps in the numerical 
simulations can vary significantly. This time mismatch 
has to be taken into consideration in order to relate the 
results of the calculations and experiments. 

The temporal variation of the applied magnetic field 
has its own time-scale, which has to be synchronized 
with the "inner clock" of the numerical simulations. This 
synchronization can be done by calculating the "real" 
time interval between two successive depinning events. 
The extremal models of the low-temperature creep are 
based on the assumption that the depinning probability 
Pd strongly depends on the maximal current density j m 
since the higher is j m the higher is the depinning proba- 
bility. We assume that this dependence has the "natural" 
exponential form 



P d oc exp 



Jm Jc 




Jm 


oc exp 


. h _ 




.31 . 



for jm < jc- (25) 



jo [t + St) = j (t ) + h Stack exp 



jm(t) 

ji 



(27) 



Changing the time scale, so that h dt tic ^ — > h we rewrite 
Eq. (|27|) in the form convenient for recursive calculations 



Jo 



jo + h exp 



Jm 
jl 



(28) 
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FIG. 7: Analogue of the current-voltage characteristics. 
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We treat now an analog of the j-E curve for the low- 
temperature magnetic flux creep. The self-organized crit- 
ically model does not introduce the electrical field E ex- 
plicitly. Therefore, we have to extend the above model 
relating the field £ to a certain characteristics of the 
model. 

When a stable critical state is established, the mag- 
netic flux redistribute inside the sample keeping the crit- 
ical current density j c almost constant. According to the 
Faraday law a varying magnetic field generates an electri- 
cal field E. The higher is the ramp rate h of the magnetic 
field, the higher is an electrical field. We assume here that 
the dependence between E and h is linear, i.e., E oc h. 
We demonstrate in Fig. how an asymptotic current 
density j is established for several values of the magnetic 
field ramp rate h. In Fig. 0is shown the dependence of 
E oc h on the asymptotic value of j, which is an analog 
of a j-E curve of Bean's critical state. This logarith- 
mic j-E curve is consistent with Eq. J3J) and numerous 
experimental data—. 



VI. SUMMARY 

We demonstrate that an extremal Robin Hood type 
model evolves to a Been's type critical state. The distri- 
bution function of the current density G(j) in this self- 
organized state was obtained by numerical simulations as 
well as analytically. We found that G(j) has a character- 
istic cut off at the critical current density. We map the 
low-temperature magnetic flux creep process to dynam- 
ics of an extremal model with Been's type critical state 
to treat magnetic flux penetration into a superconductor 
and derive an analog of the current-voltage characteris- 
tics in the flux creep regime. 
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